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Abstract 

Future climate change impacts depend on temperatures not only through changes in their 
means but also through changes in their variability. General circulation models (GCMs) pre¬ 
dict changes in both means and variability; however, GCM output should not be used directly 
as simulations for impacts assessments because GCMs do not fully reproduce present-day 
temperature distributions. This paper addresses an ensuing need for simulations of future 
temperatures that combine both the observational record and GCM projections of changes in 
means and temporal covariances. Our perspective is that such simulations should be based 
on transforming observations to account for GCM projected changes, in contrast to methods 
that transform GCM output to account for discrepancies with observations. Our method¬ 
ology is designed for simulating transient (non-stationary) climates, which are evolving in 
response to changes in CO 2 concentrations (as is the Earth at present). This work builds 
on previously described methods for simulating equilibrium (stationary) climates. Since the 
proposed simulation relies on GCM projected changes in covariance, we describe a statistical 
model for the evolution of temporal covariances in a GCM under future forcing scenarios, 
and apply this model to an ensemble of runs from one GCM, CCSM3. We find that, at least 
in CCSM3, changes in the local covariance structure can be explained as a function of the 
regional mean change in temperature and the rate of change of warming. This feature means 
that the statistical model can be used to emulate the evolving covariance structure of GCM 
temperatures under scenarios for which the GCM has not been run. When combined with an 
emulator for mean temperature, our methodology can simulate evolving temperatures under 
such scenarios, in a way that accounts for projections of changes while still retaining fidelity 
with the observational record. The emulator for variability changes is also of interest on its 
own as a summary of GCM projections of variability changes. 


1 Introduction 

Assessing the potential impacts of future climate change on areas of societal interest, such as 
agriculture and public heath, requires an understanding of how climate features important to 
those areas are expected to change. Impacts often depend on more than just the response of global 
or even local mean temperatures to greenhouse gas forcing. Many agricultural crops, for example, 
are highly sensitive to even brief periods of stress temperatures, particularly at certain times of 
the growing cycle, so crop yields can be strongly affected by changes in temperature variability 
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even in the absence of a change in mean (e.g., Wheeler et al. (2000)). In part because of examples 


like this, the climate and impacts communities have been interested in understanding changes in 
temperature variability in future climates. 

Potential future changes in temperature variability are not yet well understood. By its third 
assessment report, the Intergovernmental Panel on Climate Change (IPCC) stated that there was 
some empirical evidence for a decrease in variability at intra-annual timescales, but sparse evidence 
for changes in inter-annual variability (IPCC, 2001). More recent studies have not produced more 
definitive conclusions, with results apparently depending on specific definitions of variability and 


timescale as well as on the region being studied (IPCC, 2007). (The most recent report, IPCC 


(2013), frames variability changes in the context of extreme events, which are not a subject of 
this work.) One example of a physical mechanism that might explain variability changes at intra¬ 
annual timescales in a particular region is that changes in the polar jet stream can produce more 
persistent weather patterns over, for example, North America (e.g., Francis and Vavrus (2012)), but 


the mechanism and even detection of these changes remains controversial (Screen and Simmonds 


2013, Barnes, 2013). Implicit in the broader discussion about variability is that variability changes 


can differ by timescale of variation or geographic location. Because specific impacts will depend 
on timescale and geographic location, methods for assessing changes should be able to make such 
distinctions. That is, understanding projected changes in variability relevant to impacts is a 
problem of understanding the changes in covariance structure of a spatial-temporal held that is 
evolving in time. 

Beyond empirical studies, the primary tools used to understand and project changes in the 
distribution of climate variables are atmosphere-ocean general circulation models (GCMs). GCMs 
are deterministic, physical models that are used to generate runs of modeled climate under, for 
example, varying forcing scenarios. While GCMs are deterministic, the climate system being 
modeled is chaotic and so GCM realizations under the same forcing scenario but with different 
initial conditions will behave as if they were statistically independent. Summarizing the statistical 
properties of GCM predictions under different forcing scenarios is a challenge on its own. 

That said, GCM runs alone are often not sufficient as inputs for impacts assessments, which 
may require realistic simulations from the full distribution of the Earth’s temperatures. It is well 
understood that GCMs somewhat misrepresent observed temperature distributions under present 
conditions: regional mean temperatures may differ from observations by several degrees, and stud¬ 
ies have noted discrepancies between higher order moments of the modeled and observed climate 
distributions (IPCC (2013) and references therein). On the other hand, GCMs are assumed to pro¬ 
duce informative projections of, for example, future changes in mean temperatures due to changes 
in greenhouse gas forcing: the underlying physics are relatively realistic, and GCMs are able to 
reproduce observed temperature trends in historical forcing runs (e.g., IPCC (2013)). Projections 


of variability changes also have many consistent features across different GCMs, although current 
studies do not address changes in full covariance structures (e.g., Schneider et al. (2015); Holmes 


et al. (2015)). Impacts assessments researchers have therefore recognized a need to understand 


not only how temperature distributions are changing in GCM runs, but also how to combine those 
projections with the observational record to produce simulations of temperatures that more likely 
follow the distribution of real future temperatures. 

There are two popular classes of approaches for generating such simulations: those that modify 
GCM output to account for model-observation discrepancies (model-driven procedures), and those 
that modify observational data to account for changes projected by GCMs (observation-driven 
procedures); see, for example, Ho et al. (2012) and Hawkins et al. (2013) for reviews of common 
strategies. The most basic model-driven procedure is simple “bias correction”, where the differ- 
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ence in mean between observed temperatures and those in historical GCM realizations is assumed 
constant over time, and the estimated bias is subtracted from future GCM runs. The most basic 
observation-driven procedure is the Delta methocQ where, by contrast, changes in mean temper¬ 
ature are estimated by comparing GCM future realizations with those under historical forcing, 
and then this trend is added to the observational data. Both approaches implicitly assume that 
GCMs correctly capture changes in mean temperature; however, they can result in temperature 
simulations that have very different higher-order characteristics. 

An appealing property of observation-driven procedures like the Delta method is that they 
preserve attributes of the observations that are not explicitly accounted for in the simulation 
procedure, a property not shared by model-driven procedures. Figure[lJ top row, provides a cartoon 
illustration of this difference between the two approaches. Here, the cartoon model predicts changes 
in mean but badly misrepresents the mean and covariance structure of the observations. In such 
a setting, simple bias correction will maintain the model’s misrepresented covariance structure, 
while the Delta method yields a more realistic simulation (see Hawkins et al. (2013) for a less 
idealized example). More complicated versions of bias correction exist that attempt to correct 
for higher-order discrepancies between models and observations. Some correct for discrepancies in 


marginal distributions (e.g., Wood et al. (2004)), while others additionally attempt to correct rank 
correlation structures and inter-variable dependence structures (e.g., Piani and Haerter (2012); 


Vrac and Friederichs ( |2015 )). While such methods are more sophisticated than the simple bias 
correction illustrated in Figure [lj they too will leave intact discrepancies between the model and 
observations not accounted for in the correction procedure. If impacts assessments require realistic 
simulations from the joint distribution of temperatures across space and time, our perspective is 
that this objective is more easily met by observation-driven methods. 

Other routinely used simulation methods exist that do not fall as neatly within the model- 
driven/observation-driven dichotomy. For example, in simulations produced by stochastic weather 


generators (Semenov and Barrow, 1997, Wilks and Wilby 1999), the observations are replaced 


with synthetic data drawn from a stochastic model meant to mimic the distribution of the obser¬ 
vations. The stochastic weather generator can then be modified to account for changes predicted 
by a climate model. The drawback of this approach is that a statistical model for the observations 
is required in addition to a model for GCM projected changes, whereas observation-driven methods 
only require the latter. Synthetically generated observations will be less realistic than the observa¬ 
tions themselves. Other related methods in the statistics literature also attempt to use statistical 
models to blend information from observations with climate models to produce future simulations 
(e.g., Salazar et al. (2011)), but the proposed statistical models make very strong assumptions on 
the spatiotemporal distribution of the observations and the climate model realizations. Especially 
when projected changes from the historical climate are not very large, simulation methods should 
preserve features of the observed climate where possible. We therefore view observation-driven 
methods like the Delta method as likely to produce more realistic simulations than these methods 
as well. 

Current observational data products are not themselves without some uncertainties that will 
affect all of the simulation methods we have described, since they all rely on observational data 
in some way. For example, uncertainties emerge from the inhomogeneous distribution of observa¬ 
tion locations, measurement error of a potentially time-varying nature, and interpolation schemes 
used both to produce gridded data products and to downscale or upscale from one spatial grid to 


1 The meaning of the term “Delta method” in the geosciences, and in this work, is distinct from its typical 
use in the statistics literature for methods that employ Taylor expansions to derive asymptotic approximations of 
properties of functions of random variables. 
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Figure 1: Cartoon illustration comparing strategies for simulating temperatures that combine 
information from a model and the observational record. Columns compare simple bias correction 
(left), the Delta method (center), and our proposed method (right). Top row, the model predicts 
changes in mean temperature but no changes in variability; in this case, our proposed method 
is equivalent to the Delta method. Bottom row, the model predicts changes in both mean and 
covariance. Simple bias correction does not retain the higher order properties of the observations, 
whereas the Delta method does not account for model changes in covariance; our proposed method 
does both. 


another. The last point highlights the important issue of spatial change of support and misalign¬ 
ment (Gotway and Young, 2002): some raw observational data are point-referenced (e.g., station 
data), whereas others represent area averages (e.g., satellite data), so care must be taken to ap¬ 
propriately combine these data to produce gridded products. Some observational data products 
are distributed as ensembles in an attempt to explicitly account for uncertainties 


(e.g., 


Morice 


et al. (2012)). 1PCC (2013), Chapter 2 (Box 2.1 and elsewhere), and references therein contain 
discussions of uncertainty in the observational record. In this work, we will illustrate our method 
using one observational product, but as with all existing simulation methods that use observations, 
the method we will propose can be used with any and all available data products. 

While we have advocated for observation-driven simulation methods, an important limitation 
of the observation-driven Delta method described above, is that it does not account for changes 
in variability (Figure [lj bottom row). Some extensions of the Delta method account for changes 
in marginal variance projected by a GCM (e.g., see again Hawkins et al. (2013)), but, again, since 
variability changes need not be uniform across all timescales of variation, changes in marginal 
variance are not a complete summary of GCM projected changes in variability. Leeds et ah 


(2015) introduced an extension of the Delta method that does account for changes in the full 
temporal covariance structure projected by a GCM, but their approach is applicable only for 
equilibrated climates, in which temperatures (after preprocessing for seasonality) can be assumed 
to be stationary in time. The Earth’s climate, however, is and will continue to be in a transient 
state, in which temperatures will by definition be nonstationary in time. There is therefore an 
outstanding need for methods both to characterize changes in covariance in transient, nonstationary 
climates, and to simulate temperatures in such climates. 
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In this work, we build on the work described in Leeds et ah (2015) to develop a methodology 
for generating observation-driven simulations of temperatures in future, transient climates that 
account for transient changes in both means and temporal covariances. In Figure [lj bottom row, 
our proposed method, unlike simple bias correction or the Delta method, both accounts for the 
relevant changes projected by the cartoon model and retains other distributional properties of the 
observations. Our method reduces to the Delta method in the case that the model predicts no 
changes in variability, and reduces to the method in Leeds et al. (2015) if the past and future 
climates are both in equilibrium. Since such a simulation uses projected changes in covariances 
from a GCM, our methodology must provide a way of modeling and estimating these changes in 
transient GCM runs. The transient, nonstationary setting adds substantial challenges, and so the 
statistical modeling of changes in covariance in transient GCM runs is a primary focus of this 
paper. 

As a final complication, since GCMs are extremely computationally intensive, it is not possible 
to run a GCM under every scenario relevant for impacts assessments. In the absence of a run 
for a scenario of interest, impacts modelers may instead rely on a GCM emulator, a simpler 
procedure that produces, for example, mean temperatures that mimic what the GCM would have 
produced had it been run. Our framework for simulations can use emulated rather than true GCM 
projections. For methods that emulate mean temperatures over forcing scenarios, see |Castruccio 


et al. (2014) and references therein, and the literature stemming from the pattern scaling method 


of Santer et al. (1990). Much of the literature on climate model emulation has focused not on 


emulating model output across different forcing scenarios, but rather on emulating output with 
differing values of key climate model parameters (often for the purpose of selecting values of 


those parameters) (e.g., 

Chang et al. 

(2014 

); 

and Forest (2009); Sham Bhat et al. 

(2012 

); 


Rougier et al. 

(5 

>00< 

3); 

Sanso et al. (2008 

); 

Sanso 

Williamson et al 

• (5 

>013 

) and others). While the 


statistical concerns related to emulating climate models in parameter space are somewhat different 
from those of emulating in scenario space, the general idea remains the same: that one may use 
available climate model runs to infer properties of a run that has not been produced. In our case, 
we require an emulator for the GCM changes in covariance in addition to a mean emulator. Our 
proposed statistical model can be used for this purpose, allowing our observation-driven procedure 
to simulate future temperatures in a potentially wide range of forcing scenarios. 

The remainder of this article is organized as follows. In Section 2, we motivate and describe 
our procedure for observation-driven simulations of future temperatures in transient climates that 
accounts for projected changes in both means and temporal covariances. In Section 3, we describe 
a GCM ensemble that we use to illustrate our methodology. In Section 4, we describe a statistical 
model for the changes in temporal covariances observed in this GCM ensemble that can be used as 
an emulator for these changes; we also discuss the estimation of this statistical model. In Section 
5, we discuss results, both in terms of the GCM projections and the corresponding simulations, as 
well as the quality of our model in emulating the GCM projections. In Section 6, we give some 
concluding remarks and highlight areas for future research. 


2 Observation-driven simulations of temperatures in future 
transient climates 

Our goal is to provide a simulation of future temperatures in a transient climate under a known 
forcing scenario. In light of the preceding discussion, this simulation should reflect knowledge of 
the changes in the mean and covariance structure of future temperatures under that scenario, but 
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should otherwise preserve properties of the observed temperature record. 

Our proposed procedure is motivated by an idealization of the problem, supposing that the 
future changes in mean and temporal covariance structure are known. Following this motivation, 
we describe some modifications to the proposed procedure that we argue make the procedure more 
useful in practical settings when changes in mean and covariance must be estimated from, for 
example, GCM runs. 


2.1 Idealization 

Consider a family of multivariate (i.e., spatially referenced) Gaussian time series, z[ s \t ) at times 
t = ...,—1,0,1,... and locations l = 1,...,L, indexed by s 6 5, some set of scenarios. Write 
for the unknown mean of z[ s \t ) and assume that at each location, zf\t ) has an unknown 
evolutionary spectrum, a[ s \t,u); for details on processes with evolutionary spectra, see Priestley 
(1981). Processes with temporally varying covariance structures in general have been discussed 


extensively in the literature. An overview of a theoretical framework for understanding locally 
stationary processes, closely tied to the Priestley model, can be found in Dahlhaus (2012) and 
references therein. Our focus in this paper is on spectral methods because we view evolutionary 
spectra to be an intuitive way to characterize time-varying covariances and because the process’s 
corresponding spectral representation has useful implications for our simulation procedure. Most 
importantly, Zi(t) has, at each location and for each s, the spectral representation 


z\ S) (*) = / 4 S) (t) + [ e zult \Jal s \t,uj) d^ s) (u) , 


where ^ s \uj) is a mean zero process with orthogonal increments and unit variance; that is, 
E[clt;i 6 \uj)d{;l s \uj)*} = doj and E[d^ s \co)d^ s \uj')*] = 0 for u ^ a/, with * denoting the complex 
conjugate. Here and throughout this paper, we restrict our attention to nonstationary processes 

with evolutionary spectral representations whose transfer functions, \Ja\ s \ are real and positive. 

Suppose that we observe the time series under one scenario, z\° (f), for times t = 1,...,N 0 . 
Call Zi°\t) the observed time series; given the observed time series, we would like to generate a 
simulation of the same length as the observations, but approximately equal in distribution to an 
unobserved time series, Z/ S \t), for a given s. Since both zf\t) and z\ s \t ) are Gaussian, there is a 
class of affine transformations of z[°\t) that is equal in distribution to zf' 1 (t)\ indeed, writing Sj°' ) 
and E| ; for the covariance matrices of the observed and unobserved time series at location l, we 
have that (marginally, at each location) 




z = 


v Ah 


b + (e!”) 1/2 (s! 0) ) _1/2 ( 


,(o) 


Vi 


( 0 )> 


for any matrix square root, where x denotes the vector with entries x(t). While it is not immediately 
obvious that this fact is helpful, since the means and covariances of the two time series are unknown 
and at least for the unobserved time series cannot be directly estimated, we will describe a setting 
in which it is possible to compute this transformation (approximately) without fully knowing the 
means and covariances of the two time series. 
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The covariances of the observed and unobserved time series may be written as 
(E^) tit / = f ,u)* duj, and 


(S 


(s)\ 

l )t,t' 


j( ' t * \Ja\ s \t, uj)a\ s \t’, uj)* du, 


which follows immediately from the processes’ spectral representations. Guinness and Stein (2013) 


showed that, under some regularity conditions on the evolutionary spectra, these matrices can be 
approximated as 


do) 


-i(s) 


Cn 0 

Cn 0 


a, C 


'N 0 


-d) c [J>) 

a i I t'No I V ai 


H 

11 i , and 

H 


where H denotes the conjugate transpose and generically for some function A(t,uj) in time and 
frequency, Cjy(A) is the N x N matrix with entries 


2vr 



Cn(A),„ = ^A(t,2n(j - 

for t,j = 1,..., N. (In the setting where A is constant in time, Cn(A) is the inverse discrete time 
Fourier transform scaled by A and the result is the well-known result that the discrete time Fourier 
transform approximately diagonalizes the covariance matrix for a stationary time series observed 
at evenly-spaced intervals.) The following transformation of z[°\t) is therefore, marginally at each 
location, /, approximately equal in distribution to z\ s \t)\ 


= / 4 S) + a 


N 0 


c 


No 


.(0) 


-1 


,( 0 ) 


d 0> ). 


( 1 ) 


The crucial observation, however, is that (|Tj) can be computed exactly without fully knowing 
the means and covariances of the observed and unobserved time series. Indeed, suppose that what 
we are given are not the means and evolutionary spectra of the processes themselves, but some 
other substitute set of functions and a[ s \t,u)) satisfying, for each scenario s and at each 

location l, 

d S \t) ~ = ci, and a \ s }* ,U \ = Ww), (2) 

a) \t,u) 

for some unknown constant q and some unknown function /q(cj) that is constant in time. This 
situation is analogous to our actual predicament, where GCM runs are assumed to be more infor¬ 
mative about changes than absolute levels; indeed, one consequence of the assumptions ([2]) is that 
the substitute means and evolutionary spectra change in the same way as their true counterparts, 
so, for instance, we may write 


A ! s,0) (f) = /4 s) W - /4 u; (*) = #r(*) - and 


.( 0 ) 


(*)/ 


;( 0 ), 


p\ afi) (t,u) = 


_ a[ s \t,U)) a\ s \t,u ;) 


u) aj uj (f, oS) 


( 0 )/ 
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for the known changes in means and covariance structures. The assumption that GCM mean 
temperatures are off by a constant compared to real temperatures is essentially the assumption 
underlying both simple bias correction and the Delta method as described in Section 1, and 
we view the assumption on the evolutionary spectra as a natural extension to covariances; all 
existing simulation methods that we are aware of implicitly or explicitly make the same or similar 
assumptions (except the simple Delta method, which assumes no changes in variability at all). 

Under these assumptions, ([TJ) may be rewritten as 



A4°) + A) s ’ u; + C 


(s,0) 


N 0 





(3) 


where, to reiterate, ([TJ) and ([3]) are equal under the assumptions ([2]) because C 
C N o ( y J diag(l /hi), where diag(l jk\) is the diagonal matrix with entries l/ki(uj), so 




= C 



In light of (J3| , our proposed simulation can be computed as long as one knows (or, more realistically, 
can estimate) just the mean of the observed time series as well as the substitute evolutionary spectra 
and changes in mean. The procedure described by ([3]) is what is illustrated in the bottom right 
panel of Figure [l] 

In the case that there are no changes in covariance structure, so a\ s \t,u ) = and the 

same for their substitute counterparts, the simulation procedure ([3]) is equivalent to the Delta 
method as described in Section 1. In the case that both a\ s \t, u) and a\°\t,u) are constant 
in time, so the de-meaned time series are stationary, the procedure is the same as that described 


m 


Leeds et al. (2015). Our proposal is therefore a generalization of those two procedures, describing 


an observation-driven simulation that transforms one observed time series (possibly itself from a 
transient climate) to a simulation under a new (future, transient) scenario. 


2.2 Practical modifications to idealized procedure 


In practice, we do not actually know even substitute versions of the future changes in mean and 
covariance structure, so the procedure we have described in the preceding section must be modified 
to be made useful. 

A key assumption underlying our methodology is that GCM runs are informative about at least 
some aspect of the changes in mean and covariance structure of the real temperatures; however, 
it need not be true that the assumptions ([2]) will be satisfied by taking and a[ s \t,UJ ) to 

be the means and evolutionary spectra corresponding to real temperatures and taking fif' 1 ( t ) and 


a[ s \t,uj) to be those corresponding to temperatures under GCM runs. One possible objection to 
these assumptions is that both the observations and the GCM runs will exhibit nonstationarity in 
mean and variance due to seasonality, and it is at least plausible that the GCM representation of 
these seasonal cycles will differ from that of the observations. Leeds et al. (2015) argued that the 


seasonality can be reasonably represented as a uniformly modulated process (see Priestley (1981)) 


plus a mean seasonal cycle. That is, writing T[ s \t) for the true temperatures at time t and location 
l in scenario s, and d for the day of the year, we assume that 


T i S \t) = ^ S \t) + m T\ d ) + - A S \t)) 




















where m[ T \d) and D^ r \d) represent seasonal cycles in means and marginal variances, and z[ s \t ) 
lias mean p[ s \t) and evolutionary spectrum as above; assume a similar form for the GCM 

runs. (In the following, we will allow p\ s \t) to reflect changes in the mean seasonal cycle from 
m[ T \d) but for simplicity will assume that a^ s \t,u) has no seasonal structure; see Section 3.1 for 
details.) We estimate the seasonal cycles in mean and variability in the observations according to 


the methods described in Leeds et al. (2015); the mean seasonal cycle is modeled with the first 


ten annual harmonics and estimated via least squares, whereas the annual cycle in variability is 
estimated by a normalized moving average of windowed variances that has been averaged across 
years. We will assume that ([2]) is reasonable taking p[ s \t), a^ s \t,oj), p\ s \t), and a^(t, o>) to be the 
means and evolutionary spectra of the de-seasonalized components of real and GCM temperatures. 

Even assuming that (J2]) holds for the de-seasonalized components of observed and GCM tem¬ 
peratures, what we are given are the real and GCM temperatures themselves, not their means and 
evolutionary spectra. As such, the quantities necessary to compute ([3]) must be estimated using 
the available data. While it would be possible to estimate the corresponding evolutionary spectra 
from the GCM runs directly, note that (|3| works for any substitute function in time and frequency 
that, for each frequency, is proportional to the true evolutionary spectra at all times. If such a 
function must be estimated, there is presumably a statistical advantage to estimating a function 
that is relatively flat across frequencies. In GCM experiments, it is fairly typical to have a con¬ 
trol run under an equilibrated (often preindustrial) climate, in which at least the de-seasonalized 
temperatures can be viewed as a stationary process. Writing s = B for this equilibrated scenario, 
and a\ (uj) = d) for the corresponding spectral density of the de-seasonalized component 

of the equilibrated GCM temperatures, then if (J2]) holds, it will also be true that we can write 


Pi 




d B )(co) = d( B )(u) 


in which case cS s \t,uo) / pf' (t,uo) = a( B \u:) and i® still holds if one replaces a^ s \t,uj) with 
p\ s ' B \t,u). Moreover, we expect that the functions p\ s ' B \t,ui) will be much flatter than the 
functions af\t, c o) over the range of scenarios considered reasonable, so we expect that there should 
be some advantage in estimating these ratios rather than the evolutionary spectra themselves from 
the GCM runs. 

Writing T[ s,0 \t) for our simulation of the true temperatures under scenario s, our proposed 
procedure is therefore 


i(s,0) 


-(0) . ~(T) 

Mz + m) 


+ A i 

+diag(D; T) )C. 


(5,0) 


(4) 


'N 0 


(s,B) 


P 


f’- 1 

°iVo 


a/ p(°’ B ) 


do) 


M, 


(ob 


where x generically represents an estimate of the quantity x. In words, the procedure is, in order, 
(i) estimate and remove seasonality and mean trend in the observational record; (ii) estimate the 
future changes in mean and marginal spectra using an ensemble of GCM runs; (iii) de-correlate 
the de-trended and de-seasonalized observations using the estimated substitute function, 
obtained in step (ii); (iv) apply the estimated changes in spectra, pi s,B \ to the decorrelated data 
and invert the transformation in (iii); and (v) replace the seasonal cycles in means and variability, 
and add the new changes in mean. 

The practicality of Q depends both on our ability to obtain good estimates of all of the 
involved quantities and on our ability to compute the simulation efficiently. Gastruccio et al. (2014) 
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Figure 2: GCM ensemble [CO 2 ] trajectories. The historical scenario begins in 1870 and follows 
observed [CO 2 ] until 2010, after which it branches into the three future scenarios increasing at 
different rates until 2100. The preindustrial run maintains 1870-level [CO 2 ] until the year 4600, 
but we use only the last 1,500 years of that run. The ensemble includes eight realizations each 
under the historical, high, and low scenarios, five under the medium scenario, and one under 
preindustrial [C0 2 ]. 


and Leeds et al. (2015) collectively describe methods that essentially can be used to estimate all 
of the necessary quantities in the procedure except, crucially, the changes in evolutionary spectra. 
In the following, we will discuss modeling and estimating these functions from an ensemble of 
GCM runs. The particular statistical model we develop allows for efficient computation of the 
simulation. 


3 Description of GCM ensemble 


We study changes in the distribution of daily temperatures in an ensemble of GCM runs made with 


the Community Climate System Model Version 3 (CCSM3) (Collins et al. 2006, Yeager et al. 2006) 


at T31 atmospheric resolution (a 48 x 96 grid with a resolution of approximately 3.75° x 3.75°), 
and a 3° resolution for oceans. All runs require a long spin up; the realizations in our ensemble 
are initialized successively at ten-year intervals of the NCAR b30.048 preindustrial control run. 
Each transient realization is then forced by historical CO 2 concentrations (herein, [CO 2 ]) from 
years 1870-2010, at which point the ensemble branches into three future increasing [CO 2 ] scenarios 
for the years 2010-2100, which we name the “high”, “medium”, and “low” concentration scenarios 
(Figure |2]). For each scenario, we have a modest number of realizations (eight realizations from 
the historical, high, and low scenarios, and five realizations from the medium scenario), so the 
transient ensemble consists of about 1.1 million observations at each grid cell, or about 5 billion 
observations in total. As is typical, we will assume throughout that the ensemble members can be 
viewed as statistically independent realizations due to the system’s sensitivity to initial conditions. 

The focus of our investigation is on changes in temporal covariance structure in transient 
(nonstationary) runs of the GCM, but for the reasons described in Section 2.2, it is helpful to 
have a representation of the model’s climate in a baseline, equilibrated state. For this purpose, 
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we use a single, long run forced under preindustrial [CO 2 ] (289 ppm) for an additional 2,800 years 
past the control run initialization to ensure that the run is fully equilibrated, from which we use 
the last 1,500 years, or about 0.5 million days, for a total of about 2.5 billion observations under 
pre-industrial [C0 2 ]. 

In the following, we will index the members of the GCM ensemble by their [C0 2 ] scenario, 
s G {B,0, H, M, L}, denoting, respectively, the baseline preindustrial, historical, high, medium, 
and low scenarios. 


3.1 Data preprocessing 

The primary inferential aim of this work is to obtain estimates of p\ s ' B \t,co), the changes in 
marginal evolutionary spectra of the de-seasonalized component of daily temperatures in the GCM 
under scenario s compared to the preindustrial climate. The GCM runs have, accordingly, been 
preprocessed to remove means and seasonal cycles of variability. 

As with the observed temperatures in Section 2, we represent temperatures in the preindustrial 
run at each grid cell as a uniformly modulated process plus a mean seasonal cycle and retain 
the stationary component of this process. That is, write y\ B \t) for the temperature in the raw, 
equilibrated preindustrial GCM run at time t and location /, and again write d for the day of the 
year (the GCM does not have leap years, so d e {1,..., 365}). We represent these as 

y[ B \t) = + Di(d)x\ B \t ), 


where mi(d ) and Di(d ) are the estimated seasonal cycles in mean and marginal variance, and 
x) (t) is assumed to be stationary in time. The mean seasonal cycle and the seasonal cycle of 
marginal variance are estimated as described in Leeds et al. (2015), as also in the preceding section 
for the observational data. We retain x\ B \t ), the de-seasonalized component. 

Temperatures in the transient runs, on the other hand, will in general have an evolving mean 
in addition to nonstationarity due to seasonality. Write y\ s }(t ) for the temperature in the r’th 
realization of the transient scenario s at time t and location l, and assume the representation 

Visit) = 7 M d ) + h! A) 00 + Di(d)xft(t), 


where y\ s \t) represents an estimate of the evolving mean under scenario s (possibly including 
changes in the mean seasonal cycle) and x\ s ^ (t) is assumed to be some mean zero, but nonstationary, 
process. For simplicity we assume that the seasonal cycles of marginal variability do not evolve in 
time (overall marginal variability is still allowed to change in x\^ (t ), but such changes are assumed 
the same for each season). While the change in mean is needed for the simulation Q (see Seciton 
Al for details on its estimation), we would like to work with mean zero processes to estimate 
the changes in covariance structure. We therefore first remove from each transient realization the 
scenario average: 




1 

R, 


tvlSm, 


k =1 


where R s is the number of realizations in the ensemble under scenario s. The resulting contrasts, 
y[ s }(d), have mean zero, but still exhibit seasonal cycles in marginal variance. We thus retain the 
de-seasonalized contrasts 



Visit) 
Di(d) ' 


(5) 
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While we view each run, y\ s } (t), as independent, the contrasts, qf^ (t), are of course not independent 
across realizations within a given scenario. 

We assume that the de-seasonalized component of the preindustrial run, x[ B \t), has unknown 
marginal spectral density and the de-seasonalized component of the transient runs, x^(t), 

has unknown evolutionary spectrum a\ s \t, oS). While x\ B \t) and x fl (t) will not be equal in 
distribution to the (de-seasonalized) real-world temperatures, past or future, we assume that the 
true changes in evolutionary spectra under scenario s are equal to those of the GCM, so the 
de-seasonalized GCM and observed temperatures satisfy ([2]) and, in particular, 


(■ s,B) 

Pi 


(t,u) 





cu 


d { i s \t,ix) 




UJ 


In the following section, we discuss modeling and estimating these changes in evolutionary spectra. 


4 GCM projected changes in temporal covariance 

We describe a methodology for modeling and estimating the changes in covariance structure in a 
GCM as a function of a [CO 2 ] scenario. Our goal is not only to describe the changes in covariance in 
scenarios within our ensemble, but also to provide an emulator for the GCM changes in covariance 
for scenarios for which we have no runs. To the extent that the model we propose describes the 
GCM changes across the range of [CO 2 ] scenarios in our ensemble, the resulting emulator may be 
expected to provide good predictions of the GCM changes, at least for scenarios in some sense 
within the range spanned by our ensemble. 


4.1 A model for GCM changes in temporal covariance 


An important insight stated in Castruccio et al. (2014) is that changes in the distribution of 


temperatures in transient GCM runs under a [CO 2 ] forcing scenario should be describable in terms 
of the past trajectory of [CO 2 ]. More specifically, writing [CO 2 ](^) for the CO 2 concentration at 
time t, the distribution of temperature at time t is determined by the function 

f(t') = [C0 2 }(t-t'), for t! > 0, 

where / does not depend on t, so one does not need a different emulator for every time t. Providing 
useful statistical emulators for changes in the distribution of temperatures in transient GCM runs 
then depends on our ability to find useful functionals of the past [CO 2 ] trajectory that help explain 
those changes. 

One potentially useful summary of the past trajectory of [CO 2 ] for a given scenario is in fact 
the change in regional mean temperature relative to the preindustrial climate. We denote this 
change as A^ ,B \t) for region S. In this work, we have subdivided the T31 grid into the same 


47 regions as in Castruccio et al. (2014), chosen to be relatively homogeneous but still large 
enough to substantially reduce inter-annual variations. We estimate A^’ B ^ (t) in each region using 
a modification of the mean emulator described there (see Section Al). 

While the changes in regional mean temperature are themselves useful summaries of the past 
trajectory of [CO 2 ], it need not be true that temporal covariance structures will be the same if 
B \t') for two scenarios s and s' at two different time points t and t '. In particular, 
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the rate of change of the evolution of regional mean temperatures (d t A^’ B \t) for scenario s ) may 
capture some additional aspect of the changing climate that is also relevant for explaining changes 
in covariances. 

We have indeed found that the following model usefully describes the changes in temporal 
covariances in scenarios in our ensemble: 

log Pi S,B \tj u ) = 5joMA!s’ b) (t) + 5 n (uj)d t A { g ,B) (t). (6) 

In the case that Si 0 (oj ) and <5n(o;) are constant functions, for example, (Jb]) describes a uniformly 
modulated process. More generally, Si 0 (uj ) and <5 ji(o;) describe the patterns of changes in variability 
across frequencies associated with changes in regional mean temperature and the rate of change of 
regional mean temperature, respectively. 

Since each 5u{oS) is not scenario-dependent, model ([6]) can be thought of as an emulator for 
the GCM changes in covariance structure. That is, given an emulator for the regional mean 
temperature changes in the scenario of interest, <§ provides a prediction of the GCM changes 
in covariance structure under that scenario. In Section 5, we will discuss how well this model 
describes changes across the scenarios in our ensemble. Note that a model like (|6]) is unlikely 
to hold generically for all [CO 2 ] scenarios. In particular, such a model would be unlikely to fully 
capture the changes in variability in scenarios where [CO 2 ] changes instantaneously; such scenarios 
are typically not considered realistic. Additionally, since the changes in covariance in (Jg]) depend 
on the [CO 2 ] trajectory through the corresponding changes in mean, this model will not fully 
capture changes in variability that depend on absolute temperatures through, for example, phase 
changes between ice and water (see Figure [5] and discussion in Section 5.2). The model will 
also not fully capture GCM behavior if that behavior involved abrupt changes in the distribution 
of temperatures even under relatively smooth forcing scenarios (nonlinear responses to forcing); 
however, the model we study does not exhibit such behavior over the range of [CO 2 ] scenarios we 
study. We have found that the model is a useful description of changes in variability in scenarios 
like those in the ensemble we use here, where [CO 2 ] changes slowly and relatively smoothly over 
time and in locations not involving changing ice margins over the course of the scenario. 


4.1.1 Estimating Su(uj) 


To estimate the functions 5i 0 (u>) and 5a (<n), we adopt the intuitive approach for likelihoods for 
processes with evolutionary spectra where the usual periodogram in the Whittle likelihood is 
replaced with local periodograms over smaller blocks of time (Dahlhaus 1997). While we view 


the local Whittle likelihood approach as most suitable in our setting, several other alternative 


methods for estimating evolutionary spectra have been proposed. Neumann and Von Sachs (1997) 


use a wavelet basis expansi on; |Ombao et al.| (2002) use smooth, localized complex exponential 
basis functions; Dahlhaus (2000) proposed another likelihood approximation that replaces the 


local periodogram in the earlier work with the so-called pre-periodogram introduced by Neumann 


and Von Sachs (1997); Guinness and Stein (2013) provided an alternative generalization of the 


Whittle likelihood that they argued, at least in the settings they studied, is more accurate than 
the approximations given by either Dahlhaus (1997) or Dahlhaus (2000). However, an advantage 


of the local Whittle likelihood approach is that in addition to being intuitive, the corresponding 
score equations are computationally easier to solve in our setting when the evolutionary spectra 
evolve very slowly in time so the local periodograms can be taken over large blocks of time. 
Computation is an especially important consideration when estimating a semi-parametric model 
such as ©• Furthermore, the results from Guinness and Stein (2013) suggest that the local Whittle 
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likelihood approach may yield point estimates that are close to optimal even when the likelihood 
approximation itself is inaccurate, whereas for example they demonstrated that the approach based 
on the pre-periodogram can give unstable estimates. 

In this work, we interpret the local Whittle likelihood approach as follows. We divide each 
contrast time series, q\ s ^ (t ), defined in (J5]) and of length N s , into blocks of length M (for simplicity 
take M to be a common factor of each N s ). In our setting, we take M — 10 years, but since 
temperature variability changes very slowly over time in the scenarios we analyze, the conclusions 
are not very sensitive to the choice of M; the results are essentially the same taking M — 5 years 
or M = 30 years, for example. Upon choosing M, then for the time block, location, realization, 
and scenario indexed by 6, /, r, and s, respectively, define the local periodogram of the contrast 
time series at frequencies Uj = 0, 2n/M, ...,27r as, 



1 

27 tM 


M 

7=1 


(7) 


It is straightforward to show that the Whittle likelihood for each time block, location, and scenario 
depends on each only through the average across realizations, 



1 

Rs 


Rs 




r =1 


Likewise, for the de-seasonalized preindustrial run, define its periodogram as 


t(B) 

1 i 


(Ujj) 


1 

2ttNb 


n b 


( 8 ) 


( 9 ) 


Ill our setting, since M < Nb-, t j is defined on a coarser frequency scale than is /,! T '' fay), 
so for the purposes of estimating changes in spectra, it may be natural to aggregate the baseline 
periodogram to the coarser scale; that is, write 


f(5) 


(k-b) 


M 



{B) 


f 2?r (j + k) \ 

V n b )■ 


( 10 ) 


An approximate likelihood under model (|6j) , marginally at each location l, may then be writ¬ 
ten as the sum of the local Whittle likelihoods of the transient runs and the Whittle likelihood 
corresponding to the aggregated periodogram of the baseline run (under the usual approximation 
that the periodogram ordinates are independent at distinct Fourier frequencies): 


CiW = -)£{(*. " l)(loga< fl) (nq) + A^S ,„( UJ ) + (11) 

s,b,j 

+ _R,/^' ) (a;j)e _ ^ losa ' (fl> ^ J 7+^i, s ’ s) ' 5 ioUd+ 9 ^6 S ’ s)<5 ;iK'))| 

-yE { lo g + iJ B) (wj)/a, {B) (a;j)} • 
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where 9 = (a^ B \ S/o, <5n) and where Aj ; s,B ' > and d t A^' B ^ correspond to the values of A ( s,B \t) and 
d t A( s,B \t) for t at the midpoint of time block b. Here the ( R s — 1) factor multiplying the log- 
determinant approximation takes into account that the contrasts, q^(t), are obtained by subtract¬ 
ing off the scenario average across realizations; see Castruccio and Stein (2013) for details. 

The estimator maximizing , say 


9* = arg max Ci ( 9 ), 
6 


( 12 ) 


will yield very rough estimates of the functions a\ B \u), Sio(u), and 5/i(u;) because no smoothness 

/ D\ 

has been enforced across frequencies. The baseline spectrum, a\ '(ca), is not of particular interest 
to us, as this function is not required for the simulation ([!]). On the other hand, maximizing ( fTT| 
is clearly inadequate for estimating the functions of interest, 8iq(lo) and <5zi(a;). 

A common approach for nonparametrically estimating the spectral density of a stationary 
process is to smooth its periodogram either by kernel methods or by penalized likelihood methods. 


For estimating ratios of spectra between two stationary processes, Leeds et al. (2015) adopted a 


penalized likelihood approach whereby the penalty enforced smoothness in the ratio. Here, we opt 
to smooth the rough estimates, £ z * 0 and 8 ^, using kernel methods; that is, for i — 0,1 write as the 
final estimate for 5u 

(13) 




where Wk, h i are weights (possibly varying with j and i) satisfying Y2k w k,j,i = 1 for each i and j. 
In practice, we use weights corresponding to a kernel with a variable bandwidth that is allowed 
to decrease at lower frequencies. The reason for the variable bandwidth is that in the GCM runs 
we have analyzed, we have observed that the log ratio of spectra are typically less smooth at 
very low frequencies compared to higher frequencies. For details on the form of the weights and 
the bandwidth selection procedure we use to choose them, see Section A2. While the penalized 


likelihood approach described in Leeds et al. (2015) may be adapted for this setting, we view 


the kernel smoothing approach as more straightforward, especially when allowing for variable 
bandwidths, and have found that the approaches yield similar estimates when the bandwidth of 
the kernel is constant. 

Approximate pointwise standard errors for each 5u(uj ), and the corresponding estimate of 
logp[ b ’ B \t,ujj) may also be computed; these are described in Section A3. Having estimated our 
model, we need to compute the observation-driven simulations. Computing (jd]) efficiently is im¬ 
portant; this is described in Section A4. 

5 Results 


We estimate the model described in Section 4 using our ensemble of CCSM3 runs described in 
Section 3. In this section, we describe insights into the climate system that our estimated model 
provides and illustrate how this information is used in our proposed simulation of temperatures, 
described in Section 2. In the simulation, we use temperatures from NCEP-DOE Climate Forecast 
System Reanalysis (Saha et al. 2010) as a surrogate for observational data. The reanalysis is 


run at T62 resolution (about a 1.875° grid), which we regrid to the T31 resolution on which the 
GCM was run using an area-conserving remapping scheme. We also investigate the success of our 
proposed model in describing the changes in covariances observed in our GCM ensemble and the 
quality of our model when used as an emulator. 
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5.1 Model changes in variability 

In CCSM3, changes in variability in evolving climates can be primarily characterized using changes 
in mean temperature, with a smaller contribution by the rate of change of warming (corresponding 
to terms <5 ro and 5u in ([6]), respectively). As a consequence, the projected patterns of changes in 
variability at a given time in a given future scenario largely correspond to the patterns observed 


m 


Leeds et al. (2015): the GCM projects decreases in short timescale variability at most locations, 


but increases in longer timescale variability in some regions, especially at lower latitudes (Figure |3j 
left). 

The differences in variability between scenarios due to different rates of warming are small 
compared to the overall projected changes in variability, but also exhibit patterns. To illustrate, 
we compare changes in variability under the low scenario at year 2100 to the corresponding changes 
under the high scenario in the year of that scenario experiencing the same regional mean temper¬ 
atures as at 2100 in the low scenario (Figure |3j right); this year varies by region, ranging from 
2037 to 2044. An analogous figure is given in Section A5 that shows the estimated changes in 
variability in each of the three scenarios at years corresponding to the same change in regional 
mean temperature. In about 75% of all locations, and especially in mid- and high-latitude ocean 
locations, the changes in variability under the high scenario are larger than under the low scenario. 
Larger changes under the high scenario than under the low scenario are an indication that vari¬ 
ability changes are projected to be larger in a transient warming climate than in an equilibrated 
climate at the same temperature. 

To illustrate how these changes in covariance structure are used in our proposed simulation, 
we simulate temperatures under the high scenario at a single grid cell in the Midwestern United 
States (Figure |4j which shows the observations in 2009-2010, our simulation 89 years in the future, 
and output from one of the GCM runs in the same timeframes). Mean temperatures warm, more 
strongly in the winter than in the summer at this location, and temporal variability decreases 
overall. More specifically, variability is projected to modestly decrease at higher frequencies and 
slightly increase at low frequencies. At low frequencies, the projected log ratios are within two 
standard errors of zero, but at high frequencies are significantly smaller than zero. The extent to 
which such changes are important will of course depend on the impact domain of interest. 

The distribution of temperatures in the GCM differs strongly from that in the observations, 
with differences evident by eye in the raw time series and corresponding marginal densities. For 
example, the GCM has a stronger seasonal cycle than the observations and simulation, and greater 
variability in the winter months. See Section A5 for an additional comparisons of the space- 
time covariance structures of the observations and the GCM runs: typically, we have found that 
temperatures in nearby grid cells are more coherent in the GCM than in the observations, and 
that the coherences do not change much between the historical period and the end of the high 


scenario (Figures A4 A7). Our simulation procedure does not change the coherence structure of 
the observations. Collectively, this forms an argument for our procedure, which preserves features 
of the observations. 


5.2 Assessing model fit and quality of emulation 

To assess different aspects of how well our model describes the changes in covariance structure in 
CCSM3 in evolving climates, we show three diagnostics. First, we address in which geographic 
locations the model performs relatively better or worse. Second, we ask how well the statistical 
model performs as an emulator for a scenario on which the model has not been trained. Finally, we 
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Figure 3: Left, estimates of changes in marginal spectra, at three frequencies, for the low scenario 
at model year 2100 compared to the preindustrial climate (i.e., p\ L ' B \t, oS) at year 2100 and the 
specified frequencies). Red indicates an increase in variability and blue a decrease in variability. 
Right, estimates of (p[ H,B \t\ u)/p[ L ’ B \t, uj)) slgDl ^ u for t' corresponding to the model time under 
the high scenario with the same change in regional mean temperature as at year 2100 in the low 
scenario, where sigrp t w is the sign of the log ratio at that location, time, and frequency. Magenta 
grid cells indicate smaller changes in variability under the high scenario at the same temperature, 
whereas the cyan grid cells indicate bigger changes in variability. (Black grid cells indicate the 
roughly 5% of locations where the two estimates differ in sign, so comparing the relative magnitude 
of changes is not meaningful). Figure A3 repeats the left column maps for all the low, medium, 
and high scenarios. 


show the extent to which the rate of change of mean warming improves the quality of emulation 
over the simpler model where changes in covariance are explained solely by changes in mean 
temperature. 

To examine in which locations the model performs best and worst, we compare the deviances 
of our model at each location (Figure [5]); recall that the deviance compares the likelihood under 
our estimated model to that under the saturated model where the value of the spectrum in each 
time block, scenario, and frequency is assigned its own parameter. The deviances are largest 
at the edge of the maximum present-day sea ice extent in the Southern Ocean, and relatively 
homogeneous elsewhere. The relatively poorer fit at ice margins is expected, since variability 
decreases substantially here as sea ice retreats, and those changes are therefore based in part on 
absolute temperatures. Any statistical model based purely on changes in temperature, rather than 
absolute temperatures, will have difficulties capturing variability changes due to phase changes 
between ice and water. This result should serve as a warning against using such methods over 
locations, scenarios, and time periods in which the response to [CO 2 ] changes is highly nonlinear. 

To address how well our statistical model is able to emulate GCM projected changes in vari¬ 
ability for scenarios in some sense within the range spanned by our ensemble, we re-estimate our 
model using (a) all but the realizations under the medium scenario, and (b) only the realizations 
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Figure 4: (a), estimates of p\ H,B \t,u ) for years 2010, 2055, and 2100, at a midwestern United 
States grid cell (gray curves at ± two standard errors on the log scale), (b), part of the corre¬ 
sponding simulation computed by transforming the observational data at this grid cell via Q; the 
simulation was computed for the whole length of the observational data, but only the last two 
years are displayed, (c), a run of the GCM in the years corresponding to the simulation, (d), the 
difference between the observations and the simulation in panel (b). The overall shift upwards in 
the simulation is due to increasing mean temperature. Most of the long timescale fluctuations in 
the bottom panel are due to changes in the seasonal cycle: at this location, temperatures are pro¬ 
jected to warm more in winter months than in summer months. The short timescale fluctuations 
on the order of 0.1 degrees in the difference are due to changes in variability, which decreases in 
the future simulation, (e), marginal densities by season (labeled by corresponding months) for the 
observations in 2008-2009, the simulation from 2098-2099, the GCM from the same years under 
the historical forcing and the high future scenario scenario. 


under the medium scenario. If the conclusions we draw from (a) match those drawn from (b), 
this is evidence that we have successfully emulated the changes in covariance structure under the 
medium scenario. We compare projected changes in marginal evolutionary spectra at year 2100 
of the medium scenario, estimated under these two schemes (Figure [6]). The estimated global 
patterns of changes in variability are quite similar under the two schemes. In an absolute sense, 
the biggest differences between the two schemes are at the lowest frequencies, but recall that since 
we have reason to believe that the ratios of spectra are less smooth at lower frequencies, we smooth 
with a smaller bandwidth at those frequencies and therefore our estimates of these changes are 
more uncertain. Globally, the differences between the estimates under schemes (a) and (b) are 
within two standard errors of zero in about 60-75% of the grid cells, depending on the frequency 
of interest. (Over land, the differences are within two standard errors in about 70-80% of the grid 
cells.) The locations where there are significant differences between the two schemes are, unsur¬ 
prisingly, often those where we have argued that the model should have trouble, such as in the 
Southern Ocean. In these locations, the emulator usually underestimates variability changes. 

As discussed above, a feature of our model is that changes in variability depend not only on 
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Figure 5: Deviances, comparing the approximate likelihood under our estimated model to that 
under the saturated model where the spectrum in each time block, scenario, and frequency has 
its own parameter. The number of observations at each location is about 1.6 million days. The 
regions showing the largest deviances are those where changes in sea ice have a strong influence on 
variability; in such locations, our model based on changes in mean temperature cannot be expected 
to be a fully adequate description of changes in variability. 


the change in regional mean temperature but also on the rate at which those changes occur. One 
might ask whether the simpler model that omits the second term (i.e., 8n (ui) = 0) is just as good 
at emulating changes in variability. Figure [7] displays the predictive log likelihood ratio comparing 
the simpler model to our proposed model; by predictive log likelihood, we mean that the models 
are estimated as emulators, excluding the medium scenario realizations, and the likelihoods were 
evaluated for the medium scenario realizations (as such, no adjustment for model complexity is 
necessary). At all but three out of the 4,608 locations, the full model has a larger predictive 
likelihood than the simpler model, which indicates that for the purposes of emulation it is useful 
to allow for a nonzero term involving the rate of change change of warming. 


6 Discussion 


In this work, we describe a method for transforming observed temperatures to produce simulations 
of future temperatures when the climate is in a transient state, based on the projected changes 
in means and temporal covariances in GCM output. We believe this approach should yield more 
realistic simulations of future climate than do either GCM runs or simulations based on modifying 
GCM runs. Any observation-driven procedure is, however, of course limited by the observational 
record: for example, as described, our procedure provides exactly one simulation of future temper¬ 


atures equal in length to the observational record. As also suggested by Leeds et al. (2015), longer 


simulations could be produced either by recycling the observations entirely or by resampling them 
to generate new pseudo-observations. 

An important feature of the procedure we describe is that our simulations preserve many 
features of the observational record not accounted for explicitly in the procedure. This paper is 
concerned with changes in only the temporal covariance structure of temperatures in transient 
climates. In the methodology described here, the simulation therefore preserves, for example, 
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Figure 6: Estimates of changes in marginal spectra, at three frequencies, for the medium scenario 
at model year 2100 compared to the preindustrial climate. Left, estimates use only the medium 
scenario realizations. Right, estimates use all but the medium scenario realizations (i.e., are es¬ 
timated as an emulator). Note that the apparent sudden change in behavior for low frequency 
variability in the Pacific Ocean (top left) is an artifact of the minimum bandwidth chosen for 
smoothing in these two adjacent regions; see Section A2 for details on bandwidth selection. Loca¬ 
tions are marked with (or “x”) when the difference between the emulator and the fitted model 
is more than two (or three) standard errors away from zero. The patterns are similar under both 
schemes, with most of the differences at locations where our model would not be expected to be a 
good description of changes in variability (e.g., at ice margins). 


the spatial coherence spectra of the observations. While projected changes in spatial coherences 
in the model we study appear to be small (Section A5), changes in spatial coherences may also 
be important for societal impacts. We leave for future research the possibility of extending the 
methodology to account for such changes. This extension would be challenging and interesting, in 
part because temperatures are nonstationary in space with abrupt, local changes due to geographic 
effects. 

Another challenging and interesting extension of our methods would be to jointly simulate 
future temperature and precipitation. Our work has focused on simulating temperatures, but po¬ 
tential changes in precipitation are also important for societal impacts. While there have been some 


model-driven proposals for jointly simulating temperature and precipitation (Piani and Haerter 


2012, Vrac and Friederichs, 2015), to our knowledge most approaches (both model- and observation- 


driven) proceed by simulating the two quantities separately. Versions of the simple Delta method 
can and are used for monthly precipitation, with the assumption that the GCM captures mul¬ 
tiplicative (rather than additive) changes in rainfall amount (see, e.g., Teutschbein and Seibert 


(2012) and references therein for a review of common precipitation simulation methods). A simple 
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Figure 7: Log likelihood ratios comparing our proposed model to the reduced model taking <5u(u;) = 
0. The models were estimated using all but the medium scenario realizations and the likelihoods 
were computed with just the medium scenario realizations, so the likelihood ratio is a comparison 
of the quality of emulation (no adjustment for model complexity is needed). In all but three 
locations on the globe, the likelihood under the full model is larger, which is an indication that for 
the purposes of emulation it is useful to include 5n (cj) in the model. 


Delta method cannot, however, capture the changes on the timescale of individual rainfall events, 
whose intensity changes differently than that of time-averaged rainfall, with projections of less 
frequent but more intense storms (e.g., Trenberth (2011)). Our approach, based on spectral meth¬ 


ods, is likely also inadequate for characterizing changes in variability of daily precipitation, because 
daily precipitation often takes the value zero. More sophisticated observation-driven simulations 
methods for precipitation remains an area of research, as does the joint simulation of temperature 
and precipitation. In the context of this paper, where changes to the correlation structure of 
temperatures are detectable but not very large, we find it unlikely that separately simulating tem¬ 
perature and precipitation (as in common practice) will result in large changes to their bivariate 
dependence structure. 

The characterization of variability changes in transient climates is itself an issue of scientific 
interest. One of our motivations for developing the methods described here is to enable studying 
and ultimately comparing different GCM projections of changes in temperature variability. Most 
publicly available GCM runs (including most runs mandated by the IPCC) describe plausible 
future climates, which are necessarily in transient states. The statistical model we develop uses the 
GCM change in regional mean temperature and its rate of change to describe the GCM projected 
changes in covariance structure; in the GCM runs we study, these factors effectively summarize 
the projected changes in covariance. While we have investigated changes in variability in only one 
GCM, at relatively coarse resolution, we hope that our methods are applicable across GCMs, and 
may aid in carrying out a comparison across different GCMs in a coherent and interpretable way. 
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A1 Estimating changes in regional and local mean tempera¬ 
ture 

The model for changes in covariance, requires an estimate of the changes in regional mean 
temperature, Ag ,B \t), for scenario s. We also need changes in local mean temperature to compute 
the simulation (|4|. We estimate these using a modification of the mean emulator described in Cas- 


truccio et al. 


(2014). Write ^g\t) for the regional mean temperature at time t under scenario s 


in region S. We assume that 
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for f>s £ [0,1), and [C0 2 ]^(£) and [C0 2 ]( B ) are the C0 2 concentrations under scenario s at time 
t and under preindustrial conditions, respectively. For the harmonic terms, we take K = 6. This 
mean emulator differs from that described in Castruccio et al. (2014) in essentially two ways. First, 


here we exclude one effect in theirs that was meant to distinguish short term and long term effects 
of changes in [C0 2 ]. For smooth, monotonic scenarios like those in our ensemble, it is difficult to 
distinguish these two effects. Second, whereas their paper was concerned with emulating annual 
average temperatures, here we are interested in daily temperatures, so we need terms that account 
for the (possibly changing) mean seasonal cycle; note that the GCM runs use a year of exactly 
365 days, hence our representation of the seasonal cycle in terms of harmonics of 1/365. In most 
regions, any changes to the mean seasonal cycle are small besides an overall increase in mean, 
although in regions with strong seasonal cycles, the mean seasonal cycle tends to be damped more 
in winter months than summer months (Figure [Alj ). 

For the change in regional mean temperature used as an input in ([6]), we define A g’ B \t) = 
Pi ,sCo s '(t), which is, under the above model, the change in regional mean temperature from the 
preindustrial climate excluding changes in the seasonal cycle. For the changes in local mean 
temperature used as an input in the simulation Q, we assume that the local means are related 
to the regional means through a regional pattern scaling relationship; that is, the change in local 
mean from the preindustrial climate is taken to be proportional to the change in regional mean 
from preindustrial climate: 

hi 

This is the approach advocated by Castruccio et al. (2014]) for grid cell-level mean emulation. 


' ,(8) (t) - A B) = M/4 a) (*) - /4? : *)• 
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Figure Al: Examples of emulated (gray) and empirical (black) mean seasonal cycles in four regions 
compared between the preindustrial climate (solid) and the climate under the high scenario in the 
year 2099 (dashed). The emulations are fitted according to (Al). The empirical estimates are 
obtained by averaging across grid cells in the region of interest and across time (for the equili¬ 
brated preindustrial run) or across realizations (for the high scenario realizations). The empirical 
estimates of the 2099 seasonal cycles (red) are noisy because they represent averages across only 
eight realizations of the regional mean temperature. 


A2 Choice of variable bandwidth for estimating 5 u ( lj ) 


Here we describe the choice of weights, Wkj,i, used in the estimates of Oj) 
let Wk,j,i correspond to a variable-bandwidth quadratic kernel, 


given by (13). We 


^ k,j,l 


\o, 


k € {— Mjj, ..., Mjj} 

otherwise 


where M ] t controls the bandwidth at frequency Uj for i = 0,1, which we allow to vary as M Jtl = 
(mij — m 0t i)hj(pi ) + m 0 ,i where and m 0 , i are such that 0 < m 0ti < and 

[1(1 + cos -«))), o < w, < I 

hj{pi ) = < 1 Pi < ujj < 1 - pi 

U(l + cos(f^-f))) l-f< Uj <l, 

so m i controls the maximum bandwidth, mo,j the minimum bandwidth, and pi the frequency at 
which the bandwidths transition between and rrii^ (see Figure |A2j ). The reason we allow the 
bandwidths to decrease in frequency in this constrained way is that we have found that the log 
ratio of spectra in our ensemble are typically less smooth at very low frequencies. The specific 
form of h is somewhat arbitrary, but appears to work well in practice. 

To implement the estimator, we need to choose the parameters mp, mo,i, and pt. This is done 
via leave-one-out cross validation; that is, for each region the parameters are chosen to minimize 
the sum of the cross-validation scores for each location, 


Rcv(mi,i,m 0 ,i,Pi) = 


~ tufa)? 

(1 - W 0:j ,i ) 2 


26 










400 


350 - 

300 

250 

200 

150 

100 


1 / 

/ 

1 / 

—p = 0 

i ,y _ 

-p = 0.1 


-- p = 0.5 


■;.P-1 


0.1 0.2 0.3 0.4 

Frequency 


0.5 


Figure A2: Examples of the function Mj : i controlling the bandwidth of the kernel estimator as 
a function of the frequency being estimated. In each of these curves m 0:l = 50 (controlling the 
minimum bandwidth) and = 350 (controlling the maximum bandwidth), and the value of p, : 
(controlling the transition between the minimum and maximum bandwidths) is noted in the figure. 
The form of is such that the bandwidth of the kernel is smaller at low frequencies. 


where b* u is the maximizer of the approximate likelihood, 0 , and 5 h is the smoothed version with 
weights Wk,j,i defined above. These parameters vary by region, but typical values of (m^, mo,i,Pi) 
when estimating 8io and 8n using all runs in our ensemble are around (350, 50, 0.4) for i = 0 
and (800,400,0.7) for « = 1 (recall that the local periodograms in (11) are taken over ten-year 
blocks, so the raw estimates of 5[ t are calculated at 1,825 frequencies). Since the parameters vary 
regionally, there may be sudden changes in estimates of the changes in spectra near boundaries of 
the region (see Figure [6j top right). One could consider a post-hoc spatial smoothing of either the 
bandwidths or the changes in spectra to avoid this, but this is not explored in this paper. 


A3 Standard errors for estimates of 5 u ( lj ) 

Under ( JIl| , the Fisher information for (<5; 0; 8n, a ^) is a block diagonal matrix with blocks for each 
frequency, uij, approximately (under standard approximations for periodograms) equal to 
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where R' s = R s — 1. We calculate standard errors under the usual approximation that the variance 
of the maximum likelihood estimate§ is the inverse information matrix; such an approximation 
can be justified by, for example, considering either the number of independent GCM runs or the 
number of time blocks to be large. The resulting covariance matrix for 5f 0 (ojj) and (uq-), the 
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Then the variance of the smoothed estimates, 5, dehned in (13), at frequencies not close to 0 or it 

is 


Vi i Vi i 

d lO, d lO WO^Zl 

V Sn,Sn - 

VsioAo Sfc w l,j ,0 Vs l0 ,s n WkjflWkj.l 

V SnAiJ2k w l,j,l 

When 0 or n is within the local bandwidth for the frequency of interest, we make the correction 
to the above to account for the fact that (ujj ) is periodic and even (here omitted for simplicity). 
The variance of log p\ s,B \t , u>), the smoothed estimate of log ratio of spectra, is then approximately 

Varlogp<‘ ,B, (f,w) ss + (d,A , S ‘' B> (t)) 2 V inAl 

+2Af B) mA l V(t)V Sl an 



A4 Computing simulations 


Computing Q efficiently requires the ability to quickly compute the products Cn t y-\/p( s ’ B ')j x 
and x for a vector x. The matrix-vector products may be written as 




e iL)jt 


m^ AuJit 


where the approximation just truncates the Taylor series at P — 1, then uses the binomial expansion 
and changes the order of summation. This approximation is the weighted sum of P(P + l)/2 
inverse discrete Fourier transforms; if P can be taken to be only modestly large so that the Taylor 
approximation is accurate, this sum can be computed efficiently. In our setting, the changes 
in variability projected by the GCM are small and P need not be very large in order for the 
approximation to be quite good; we have found that taking P = 10 is more than enough to give 
accurate approximations given the magnitude of the estimated changes in variability. 

To compute the matrix-inverse vector products, we have found, as in Guinness and Stein (2013), 
that iterative methods work well. In order to work well, these require the ability to compute forward 
multiplication quickly, which we have just described, as well as a good preconditioner; since the 
historical time series is only mildly nonstationary, we have found that the Fourier transform scaled 
by the square root of the average of p^ ,B \t^u) over time works well as a preconditioner. 


Guinness and Stein (2013 
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A5 Additional supplementary figures 



Figure A3: Estimates of changes in marginal spectra, at three frequencies, at the year correspond¬ 
ing to that with the same regional mean temperature as in the low scenario in model year 2100 
(i.e., left, p[ L ’ B \t,u}) at year 2100; middle, p\ M ’ B \t,u>) at years ranging from 2050-2056 depending 
on region; right, p\ H ' B \t,uj) at years ranging from 2037-2044 depending on region). Red indicates 
an increase in variability and blue a decrease in variability. The estimated differences between 
scenarios are small because the effect of the rate of change of warming on the estimated changes 
in variability are smaller than the effect of the overall regional mean change in temperature. A 
direct comparison of the estimates between the high and low scenario is shown in Figure [3] of the 
main text. 
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Figure A4: Estimates of coherence spectra between temperatures in the Midwestern grid cell 
discussed in Section 5.1 and its four neighboring grid cells; the coordinates corresponding to each 
neighbor are indicated (and figures are positioned geographically). Estimates are shown for the 
ten-year blocks 1980-1989, 1990-1999, and 2000-2009 as well as the corresponding time 90 years 
in the future under the high scenario (i.e., 2070-2079, 2080-2089, and 2090-2099). Estimates 
in red correspond to temperatures from CCSM3 under historical forcing (solid) and in the high 
scenario future (dashed). Those in blue correspond to temperatures from NCEP CFSR and the 
simulation (the coherences do not change in the simulation). The time series were processed to 
remove means and (marginal) seasonal variability as described in the main text. Marginal spectra, 
co-, and quadrature spectra were estimated by averaging the raw estimates for each run (for the 
CCSM3 runs, where we have multiple runs) and then smoothing each component using a quadratic 
kernel; since we cannot average over multiple runs of reanalysis, the estimates for reanalysis were 
smoothed using a slightly wider bandwidth. While the coherence spectra are arguably similar in 
the East-West direction, CCSM3 has temperatures that are more coherent than the reanalysis’s 
in the North-South direction. Coherences do not change much in the future GCM run, which 
suggests that our approach for the simulations may be adequate. 
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Figure A5: The same as Figure |A4| but for a central grid cell in Northern Europe. Temperatures 
are more coherent in CCSM3 than in the reanalysis. 
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Figure A6: The same as Figure |A4| but for a central grid cell in East Asia. Temperatures are more 
coherent in CCSM3 than in the reanalysis. 
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Figure A7: The same as Figure [A4| but for a central grid cell in the Sahara. Here temperatures 
are more coherent in CCSM3 than in reanalysis at the higher frequencies. 
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Figure A8: The locations of the central grid cells corresponding to the coherency plots in Figures |A4| 


A7 (labeled by figure number). The grid cell for Figure A4 corresponds to the one discussed in 


Section 5.1 of the main text. 
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